home *** CD-ROM | disk | FTP | other *** search
/ Games of Daze / Infomagic - Games of Daze (Summer 1995) (Disc 1 of 2).iso / x2ftp / msdos / math / newmat08 / bandmat.cpp < prev    next >
Encoding:
C/C++ Source or Header  |  1995-01-18  |  10.7 KB  |  398 lines
  1. //$$ bandmat.cpp                     Band matrix definitions
  2.  
  3. // Copyright (C) 1991,2,3,4: R B Davies
  4.  
  5. #define WANT_MATH                    // include.h will get math fns
  6.  
  7. #include "include.h"
  8.  
  9. #include "newmat.h"
  10. #include "newmatrc.h"
  11.  
  12. //#define REPORT { static ExeCounter ExeCount(__LINE__,10); ++ExeCount; }
  13.  
  14. #define REPORT {}
  15.  
  16.  
  17.  
  18.  
  19.  
  20.  
  21.  
  22.  
  23.  
  24.  
  25.  
  26.  
  27. BandMatrix::BandMatrix(const BaseMatrix& M)
  28. {
  29.    REPORT // CheckConversion(M);
  30.    MatrixConversionCheck mcc;
  31.    GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::BM);
  32.    GetMatrix(gmx); CornerClear();
  33. }
  34.  
  35. void BandMatrix::SetParameters(const GeneralMatrix* gmx)
  36. {
  37.    MatrixBandWidth bw = gmx->BandWidth();
  38.    lower = bw.lower; upper = bw.upper;
  39. }
  40.  
  41. void BandMatrix::ReDimension(int n, int lb, int ub)
  42. {
  43.    REPORT
  44.    Tracer tr("BandMatrix::ReDimension");
  45.    if (lb<0 || ub<0) Throw(ProgramException("Undefined bandwidth"));
  46.    lower = (lb<=n) ? lb : n-1; upper = (ub<=n) ? ub : n-1;
  47.    GeneralMatrix::ReDimension(n,n,n*(lower+1+upper)); CornerClear();
  48. }
  49.  
  50. void BandMatrix::operator=(const BaseMatrix& X)
  51. {
  52.    REPORT // CheckConversion(X);
  53.    MatrixConversionCheck mcc;
  54.    Eq(X,MatrixType::BM); CornerClear();
  55. }
  56.  
  57. void BandMatrix::CornerClear() const
  58. {
  59.    // set unused parts of BandMatrix to zero
  60.    REPORT
  61.    int i = lower; Real* s = store; int bw = lower + 1 + upper;
  62.    while (i)
  63.       { int j = i--; Real* sj = s; s += bw; while (j--) *sj++ = 0.0; }
  64.    i = upper; s = store + storage;
  65.    while (i)
  66.       { int j = i--; Real* sj = s; s -= bw; while (j--) *(--sj) = 0.0; }
  67. }
  68.  
  69. MatrixBandWidth MatrixBandWidth::operator+(const MatrixBandWidth& bw) const
  70. {
  71.    int l = bw.lower; int u = bw.upper;
  72.    l = (lower < 0 || l < 0) ? -1 : (lower > l) ? lower : l;
  73.    u = (upper < 0 || u < 0) ? -1 : (upper > u) ? upper : u;
  74.    return MatrixBandWidth(l,u);
  75. }
  76.  
  77. MatrixBandWidth MatrixBandWidth::operator*(const MatrixBandWidth& bw) const
  78. {
  79.    int l = bw.lower; int u = bw.upper;
  80.    l = (lower < 0 || l < 0) ? -1 : lower+l;
  81.    u = (upper < 0 || u < 0) ? -1 : upper+u;
  82.    return MatrixBandWidth(l,u);
  83. }
  84.  
  85. MatrixBandWidth MatrixBandWidth::minimum(const MatrixBandWidth& bw) const
  86. {
  87.    int l = bw.lower; int u = bw.upper;
  88.    if ((lower >= 0) && ( (l < 0) || (l > lower) )) l = lower;
  89.    if ((upper >= 0) && ( (u < 0) || (u > upper) )) u = upper;
  90.    return MatrixBandWidth(l,u);
  91. }
  92.  
  93. UpperBandMatrix::UpperBandMatrix(const BaseMatrix& M)
  94. {
  95.    REPORT // CheckConversion(M);
  96.    MatrixConversionCheck mcc;
  97.    GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::UB);
  98.    GetMatrix(gmx); CornerClear();
  99. }
  100.  
  101. void UpperBandMatrix::operator=(const BaseMatrix& X)
  102. {
  103.    REPORT // CheckConversion(X);
  104.    MatrixConversionCheck mcc;
  105.    Eq(X,MatrixType::UB); CornerClear();
  106. }
  107.  
  108. LowerBandMatrix::LowerBandMatrix(const BaseMatrix& M)
  109. {
  110.    REPORT // CheckConversion(M);
  111.    MatrixConversionCheck mcc;
  112.    GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::LB);
  113.    GetMatrix(gmx); CornerClear();
  114. }
  115.  
  116. void LowerBandMatrix::operator=(const BaseMatrix& X)
  117. {
  118.    REPORT // CheckConversion(X);
  119.    MatrixConversionCheck mcc;
  120.    Eq(X,MatrixType::LB); CornerClear();
  121. }
  122.  
  123. BandLUMatrix::BandLUMatrix(const BaseMatrix& m)
  124. {
  125.    REPORT
  126.    Tracer tr("BandLUMatrix");
  127.    GeneralMatrix* gm = ((BaseMatrix&)m).Evaluate(MatrixType::BM);
  128.    m1 = ((BandMatrix*)gm)->lower; m2 = ((BandMatrix*)gm)->upper;
  129.    GetMatrix(gm);
  130.    if (nrows!=ncols) Throw(NotSquareException(*this));
  131.    d = TRUE; sing = FALSE;
  132.    indx = new int [nrows]; MatrixErrorNoSpace(indx);
  133.    MONITOR_INT_NEW("Index (BndLUMat)",nrows,indx)
  134.    storage2 = nrows * m1;
  135.    store2 = new Real [storage2]; MatrixErrorNoSpace(store2);
  136.    MONITOR_REAL_NEW("Make (BandLUMat)",storage2,store2)
  137.    ludcmp();
  138. }
  139.  
  140. BandLUMatrix::~BandLUMatrix()
  141. {
  142.    MONITOR_INT_DELETE("Index (BndLUMat)",nrows,indx)
  143.    MONITOR_REAL_DELETE("Delete (BndLUMt)",storage2,store2)
  144. #ifdef Version21
  145.    delete [] indx; delete [] store2;
  146. #else
  147.    delete [nrows] indx; delete [storage2] store2;
  148. #endif
  149. }
  150.  
  151. MatrixType BandLUMatrix::Type() const { return MatrixType::BC; }
  152.  
  153.  
  154. LogAndSign BandLUMatrix::LogDeterminant() const
  155. {
  156.    if (sing) return 0.0;
  157.    Real* a = store; int w = m1+1+m2; LogAndSign sum; int i = nrows;
  158.    while (i--) { sum *= *a; a += w; }
  159.    if (!d) sum.ChangeSign(); return sum;
  160. }
  161.  
  162. GeneralMatrix* BandMatrix::MakeSolver()
  163. {
  164.    REPORT
  165.    GeneralMatrix* gm = new BandLUMatrix(*this);
  166.    MatrixErrorNoSpace(gm); gm->ReleaseAndDelete(); return gm;
  167. }
  168.  
  169.  
  170. void BandLUMatrix::ludcmp()
  171. {
  172.    REPORT
  173.    Real* a = store;
  174.    int i = m1; int j = m2; int k; int n = nrows; int w = m1 + 1 + m2;
  175.    while (i)
  176.    {
  177.       Real* ai = a + i;
  178.       k = ++j; while (k--) *a++ = *ai++;
  179.       k = i--; while (k--) *a++ = 0.0;
  180.    }
  181.  
  182.    a = store; int l = m1;
  183.    for (k=0; k<n; k++)
  184.    {
  185.       Real x = *a; i = k; Real* aj = a;
  186.       if (l < n) l++;
  187.       for (j=k+1; j<l; j++)
  188.          { aj += w; if (fabs(x) < fabs(*aj)) { x = *aj; i = j; } }
  189.       indx[k] = i;
  190.       if (x==0) { sing = TRUE; return; }
  191.       if (i!=k)
  192.       {
  193.          d = !d; Real* ak = a; Real* ai = store + i * w; j = w;
  194.          while (j--) { x = *ak; *ak++ = *ai; *ai++ = x; }
  195.       }
  196.       aj = a + w; Real* m = store2 + m1 * k;
  197.       for (j=k+1; j<l; j++)
  198.       {
  199.          *m++ = x = *aj / *a; i = w; Real* ak = a;
  200.      while (--i) { Real* aj1 = aj++; *aj1 = *aj - x * *(++ak); }
  201.          *aj++ = 0.0;
  202.       }
  203.       a += w;
  204.    }
  205. }
  206.  
  207. void BandLUMatrix::lubksb(Real* B, int mini)
  208. {
  209.    REPORT
  210.    Tracer tr("BandLUMatrix::lubksb");
  211.    if (sing) Throw(SingularException(*this));
  212.    int n = nrows; int l = m1; int w = m1 + 1 + m2;
  213.  
  214.    for (int k=0; k<n; k++)
  215.    {
  216.       int i = indx[k];
  217.       if (i!=k) { Real x=B[k]; B[k]=B[i]; B[i]=x; }
  218.       if (l<n) l++;
  219.       Real* m = store2 + k*m1; Real* b = B+k; Real* bi = b;
  220.       for (i=k+1; i<l; i++)  *(++bi) -= *m++ * *b;
  221.    }
  222.  
  223.    l = -m1;
  224.    for (int i = n-1; i>=mini; i--)
  225.    {
  226.       Real* b = B + i; Real* bk = b; Real x = *bk;
  227.       Real* a = store + w*i; Real y = *a;
  228.       int k = l+m1; while (k--) x -=  *(++a) * *(++bk);
  229.       *b = x / y;
  230.       if (l < m2) l++;
  231.    }
  232. }
  233.  
  234. void BandLUMatrix::Solver(MatrixRowCol& mcout, const MatrixRowCol& mcin)
  235. {
  236.    REPORT
  237.    Real* el = mcin.store; int i = mcin.skip;
  238.    while (i--) *el++ = 0.0;
  239.    el += mcin.storage; i = nrows - mcin.skip - mcin.storage;
  240.    while (i--) *el++ = 0.0;
  241.    lubksb(mcin.store, mcout.skip);
  242. }
  243.  
  244. // Do we need check for entirely zero output?
  245.  
  246.  
  247. void UpperBandMatrix::Solver(MatrixRowCol& mcout,
  248.    const MatrixRowCol& mcin)
  249. {
  250.    REPORT
  251.    Real* elx = mcin.store+mcout.skip; int i = mcin.skip-mcout.skip;
  252.    while (i-- > 0) *elx++ = 0.0;
  253.    int nr = mcin.skip+mcin.storage; elx = mcin.store+nr; Real* el = elx;
  254.    int j = mcout.skip+mcout.storage-nr; i = nr-mcout.skip;
  255.    while (j-- > 0) *elx++ = 0.0;
  256.  
  257.    Real* Ael = store + (upper+1)*(i-1)+1; j = 0;
  258.    while (i-- > 0)
  259.    {
  260.       elx = el; Real sum = 0.0; int jx = j;
  261.       while (jx--) sum += *(--Ael) * *(--elx);
  262.       elx--; *elx = (*elx - sum) / *(--Ael);
  263.       if (j<upper) Ael -= upper - (++j); else el--;
  264.    }
  265. }
  266.  
  267. void LowerBandMatrix::Solver(MatrixRowCol& mcout,
  268.    const MatrixRowCol& mcin)
  269. {
  270.    REPORT
  271.    Real* elx = mcin.store+mcout.skip; int i = mcin.skip-mcout.skip;
  272.    while (i-- > 0) *elx++ = 0.0;
  273.    int nc = mcin.skip; i = nc+mcin.storage; elx = mcin.store+i;
  274.    int nr = mcout.skip+mcout.storage; int j = nr-i; i = nr-nc;
  275.    while (j-- > 0) *elx++ = 0.0;
  276.  
  277.    Real* el = mcin.store+nc; Real* Ael = store + (lower+1)*nc + lower; j = 0;
  278.    while (i-- > 0)
  279.    {
  280.       elx = el; Real sum = 0.0; int jx = j;
  281.       while (jx--) sum += *Ael++ * *elx++;
  282.       *elx = (*elx - sum) / *Ael++;
  283.       if (j<lower) Ael += lower - (++j); else el++;
  284.    }
  285. }
  286.  
  287.  
  288. LogAndSign BandMatrix::LogDeterminant() const
  289. {
  290.    REPORT
  291.    BandLUMatrix C(*this); return C.LogDeterminant();
  292. }
  293.  
  294. LogAndSign LowerBandMatrix::LogDeterminant() const
  295. {
  296.    REPORT
  297.    int i = nrows; LogAndSign sum; Real* s = store + lower; int j = lower + 1;
  298.    while (i--) { sum *= *s; s += j; }
  299.    ((GeneralMatrix&)*this).tDelete(); return sum;
  300. }
  301.  
  302. LogAndSign UpperBandMatrix::LogDeterminant() const
  303. {
  304.    REPORT
  305.    int i = nrows; LogAndSign sum; Real* s = store; int j = upper + 1;
  306.    while (i--) { sum *= *s; s += j; }
  307.    ((GeneralMatrix&)*this).tDelete(); return sum;
  308. }
  309.  
  310. GeneralMatrix* SymmetricBandMatrix::MakeSolver()
  311. {
  312.    REPORT
  313.    GeneralMatrix* gm = new BandLUMatrix(*this);
  314.    MatrixErrorNoSpace(gm); gm->ReleaseAndDelete(); return gm;
  315. }
  316.  
  317. SymmetricBandMatrix::SymmetricBandMatrix(const BaseMatrix& M)
  318. {
  319.    REPORT  // CheckConversion(M);
  320.    MatrixConversionCheck mcc;
  321.    GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::SB);
  322.    GetMatrix(gmx);
  323. }
  324.  
  325. GeneralMatrix* SymmetricBandMatrix::Transpose(TransposedMatrix*, MatrixType mt)
  326. { REPORT  return Evaluate(mt); }
  327.  
  328. LogAndSign SymmetricBandMatrix::LogDeterminant() const
  329. {
  330.    REPORT
  331.    BandLUMatrix C(*this); return C.LogDeterminant();
  332. }
  333.  
  334. void SymmetricBandMatrix::SetParameters(const GeneralMatrix* gmx)
  335. { lower = gmx->BandWidth().lower; }
  336.  
  337. void SymmetricBandMatrix::ReDimension(int n, int lb)
  338. {
  339.    REPORT
  340.    Tracer tr("SymmetricBandMatrix::ReDimension");
  341.    if (lb<0) Throw(ProgramException("Undefined bandwidth"));
  342.    lower = (lb<=n) ? lb : n-1;
  343.    GeneralMatrix::ReDimension(n,n,n*(lower+1));
  344. }
  345.  
  346. void SymmetricBandMatrix::operator=(const BaseMatrix& X)
  347. {
  348.    REPORT // CheckConversion(X);
  349.    MatrixConversionCheck mcc;
  350.    Eq(X,MatrixType::SB);
  351. }
  352.  
  353. void SymmetricBandMatrix::CornerClear() const
  354. {
  355.    // set unused parts of BandMatrix to zero
  356.    REPORT
  357.    int i = lower; Real* s = store; int bw = lower + 1;
  358.    while (i)
  359.       { int j = i--; Real* sj = s; s += bw; while (j--) *sj++ = 0.0; }
  360. }
  361.  
  362. MatrixBandWidth SymmetricBandMatrix::BandWidth() const
  363.    { return MatrixBandWidth(lower,lower); }
  364.  
  365. inline Real square(Real x) { return x*x; }
  366.  
  367.  
  368. Real SymmetricBandMatrix::SumSquare() const
  369. {
  370.    REPORT
  371.    CornerClear();
  372.    Real sum1=0.0; Real sum2=0.0; Real* s=store; int i=nrows; int l=lower;
  373.    while (i--)
  374.       { int j = l; while (j--) sum2 += square(*s++); sum1 += square(*s++); }
  375.    ((GeneralMatrix&)*this).tDelete(); return sum1 + 2.0 * sum2;
  376. }
  377.  
  378. Real SymmetricBandMatrix::SumAbsoluteValue() const
  379. {
  380.    REPORT
  381.    CornerClear();
  382.    Real sum1=0.0; Real sum2=0.0; Real* s=store; int i=nrows; int l=lower;
  383.    while (i--)
  384.       { int j = l; while (j--) sum2 += fabs(*s++); sum1 += fabs(*s++); }
  385.    ((GeneralMatrix&)*this).tDelete(); return sum1 + 2.0 * sum2;
  386. }
  387.  
  388. Real SymmetricBandMatrix::Sum() const
  389. {
  390.    REPORT
  391.    CornerClear();
  392.    Real sum1=0.0; Real sum2=0.0; Real* s=store; int i=nrows; int l=lower;
  393.    while (i--)
  394.       { int j = l; while (j--) sum2 += *s++; sum1 += *s++; }
  395.    ((GeneralMatrix&)*this).tDelete(); return sum1 + 2.0 * sum2;
  396. }
  397.  
  398.